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We demonstrate how the physics of multiboson correlation interference leads to the computational 
complexity of linear optical interferometers based on correlation measurements in the degrees of 
freedom of the input bosons. In particular, we address the task of MultiBoson Correlation Sampling 
(MBCS) from the probability distribution associated with polarization- and time-resolved detections 
at the output of random linear optical networks. We show that the MBCS problem is fnndamentally 
hard to solve classically even for nonidentical input photons, regardless of the color of the photons, 
making it also very appealing from an experimental point of view. These results fully manifest the 
qnantum computational supremacy inherent to the fnndamental nature of quantum interference. 


Motivation. The interference of multiple bosons 
based on high-order correlation measurements [1-3] in 
a linear network is a phenomenon that is fundamental 
in atomic, molecular, and optical physics. The richness 
of its features gives rise to a wide variety of applications 
in quantum information processing [1, 4, 5], quantum 
metrology [6-8], and imaging [9[. Already correlated de¬ 
tections of two bosons after the interaction with a bal¬ 
anced beam splitter reveal an interference effect of truly 
quantum mechanical origin [10-13] : both particles always 
end up in the same output port due to the destructive in¬ 
terference of the two-boson quantum paths in which the 
bosons are either both reflected or both transmitted. 

Going to higher-order correlation measurements in op¬ 
tical networks of large dimensions, multiboson interfer¬ 
ence becomes increasingly complex, promising a compu¬ 
tational power that is not achievable classically [14, 15[. 
Multiphoton correlation experiments with more than two 
photons have already been performed [16-26], providing 
an important milestone towards experiments of higher 
orders [27, 28] . 

These experiments are usually based on joint measure¬ 
ments at the interferometer output ports “classically” av¬ 
eraging over the photons’ degrees of freedom (e.g. time, 
polarization). In this context, Aaronson and Arkhipov 
argued the computational hardness of multiboson inter¬ 
ference in linear optics for identical bosons by introduc¬ 
ing the well-known boson sampling problem [14] . Does 
this computational hardness also occur for nonidentical 
photons? While the computational complexity for par¬ 
tially distinguishable photons is still not known [15], it is 
clear that boson sampling becomes computationally triv¬ 
ial for fully distinguishable photons when the information 
about the detection times and polarizations is completely 
ignored. 

However, recent technological advances have enabled 
experimentalists to produce arbitrarily polarized single 
photons with near arbitrary spectral and temporal prop¬ 
erties [29-31] which can be “read out” by time- and 
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polarization-resolving measurements [1, 32-35] with ex- 
tremely fast detectors [36] . This makes it possible to 
encode entire “quantum alphabets” in the degrees of free¬ 
dom of multiple photons [37, 38] and to retrieve the en¬ 
coded information by correlation measurements in those 
degrees of freedom, representing a valuable tool in quan¬ 
tum information processing [1, 39-52]. 

All these remarkable technological achievements now 
allow experimentalists to fully address the following fun¬ 
damental questions about the interplay between the 
physics and the complexity of multiboson interference: 
How do the spectral distributions of N nonidentical pho¬ 
tons determine the occurrence of V-photon interference 
events in time- and polarization-resolving correlated mea¬ 
surements? How and to what degree is this occurrence 
connected with computational complexity? Does compu¬ 
tational hardness really disappear for input bosons that 
are completely distinguishable in their spectra? This let¬ 
ter aims to answer all these important questions, from 
both a fundamental and an experimental point of view, 
demonstrating the inherent computational complexity of 
the physics of multiboson correlation interference even 
for nonidentical photons. 

MultiBoson Correlation Sampling (MBCS). 

We consider N single photons prepared at the N input 
ports of a linear interferometer (see Fig. 1) with 2M ^ 
2N ports. The interferometer unitary transformation U 
is chosen randomly according to the Haar measure and 
is implemented by using a polynomial number (in M) 
of passive linear optical elements [53]. The state of N 
single photons injected in a set S of N input ports s G S 
is given by 

|5) :=(g)|l[a«(8)|0)., 

s^S 

with the single photon states 

OO 

’= X (^A •€s(w))al,A(‘*^)|0)s> (1) 

A=l,2^ 

where {61,62} is an arbitrary polarization basis and 
Uj a(^) creation operator for the frequency mode 
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FIG. 1. General setup for multiboson correlation sampling. 
N single photons are injected into an Ai-port subset S of the 
M ^ N input ports of a random linear interferometer. At 
the output of the interferometer, they are detected in one 
of the possible port samples V containing N of the M out¬ 
put ports at corresponding detection times and polarizations 
{td,Pd}deT>- For each output port sample D and given input 
configuration 5, the evolution through the interferometer is 
fully described hy a. N x N submatrix of the M x M 

interferometer matrix U. 


oj and the polarization A [54] . The complex spectral am¬ 
plitude 


|s(a;) := Vs - Ws) (2) 

is defined by the spectral shape ^s(w —Ws) S E, (centered 
around the central frequency (photon color) Ws and with 
normalization f dio |^s(a;)| = 1 ), the polarization 

and the time tos of emission of the photon injected in 
the port s S 5. For simplicity, we consider input-photon 
spectra satisfying the narrow bandwidth approximation 
and a polarization-independent interferometric evolution 
with equal propagation time At for each possible path 
from an input source to a detector at the interferometer 
output. 

Given such a multiboson interferometer and assuming 
identical photons, ^ ys £ S, the boson sampling 

problem [14] was defined by Aaronson and Arkhipov as 
the task of sampling from the probability distribution 
over the output port samples T), regardless of detection 
times and polarizations. We address here an interesting 
generalization of this famous problem by introducing the 
problem of MultiBoson Correlation Sampling (MBCS) 
[1, 55]. The MBCS problem is defined as the task of 
sampling at the interferometer output from the probabil¬ 
ity distribution associated with time- and polarization¬ 
resolving correlation measurements. Each possible sam¬ 
ple corresponds to an Af-photon detection event at an N- 
port subset T) of the M output ports at given times and 
polarizations {tdTPd\d^v, with pd G { 61 , 62 } [56[. The 
Al-photon detection probability rate corresponding to a 
sample (P, {f^, pdjdg-p) depends [1] on both the N x N 
submatrix 


of the M X M unitary matrix U describing the interfer¬ 
ometer, and the Fourier transforms 

Xs{t) ■■= At) 

= VsXs{t-tos-At)e^^^^*-*»‘-^*^ (3) 

of the single-photon spectra ^s(a;) in Eq. (2) (with Xsit) 
being the Fourier transform of ^^(a;)). Defining the ma¬ 
trices 


'^{£’2} (Pd ■ Xsitd))]dev 

and using the definition 


N 

permAd := ^ 

i—1 


of the permanent of a matrix Ad, where the sum runs over 
all permutations a in the symmetric group Ejy, the prob¬ 
ability rate of an Al-fold detection event {V, {td,Pd}d£v) 
is 


r^'D,s) 

'^{G.Pdl 


nermT*'®’*^^ 
P®™ '{td,Pd} 


( 4 ) 


for ideal photodetectors. 

By considering an integration time T/ short enough 
such that 


Vfd : Xs{t-tos-At)xs'{t-tos>-At)e^^‘^‘ « const. 

yt G [td — Tj,td + Ti],\/s, s' G S, (5) 

we obtain, for a detection sample (2?, {fd,Pd}dGX>)) the 
probability 


42,2i = (2^^)"|pe-r£2>r (6) 

of an A^-fold detection in the time intervals {[id — Ti,td + 
’7/]}dG'Dj where the detection time axes are discretized 
with step width 2Ti. 

We emphasize that, for each possible sample 
{'D,{td,Pd}dGTi), the probability in Eq. ( 6 ) is at most 
exponentially small in N, as demonstrated in Theorem 1 
in the Supplemental Material [57] . 

Exact MBCS. Obviously, the complexity of sam¬ 
pling exactly from the probability distribution defined 
by Eq. (6) depends on the Al-tuples {^s}sg 5 of single¬ 
photon input spectra in Eq. (2) [58]. 

With that in mind, in order to establish the complex¬ 
ity of exact MBCS, it is useful to define the Al-photon 
interference matrix with elements 

00 

a(s,s') J dt \xsit - tos)\\Xs'(,t - tos')\ < I, 

— 00 

( 7 ) 


:= [Ud,sWv 

sGS 


with s, s' G S, depending on the pairwise overlaps of the 
absolute values of the temporal single-photon detection 
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amplitudes [59] Xs {t — tos — ^t) and of the 

polarizations Vg in Eq. (3) . For non-vanishing elements 

0 < a(s, s') < 1 Vs, s' S S, ( 8 ) 

there exists a time interval T and at least a polarization 
G { 61 , 62 }, such that 

ex-Xs(id)^0 Vtd € T,Vs € S,Vd € V. 


It is then ensured that for each detection sample 
(V,{td,Pd}dev), with td GT,pd = e-^ \/d G V, the input 
photons are indistinguishable at the detectors: this leads 
to the interference of all possible N\ 7V-photon quantum 
paths manifested by the coherent superposition of all cor¬ 
responding, non-vanishing N\ 7V-photon detection ampli¬ 
tudes in Eq. (4). Therefore, only the conditions (5) and 
(8) for the nonidentical input spectra }sg 5 in Eq- (2) 
are enough to ensure the occurrence of 7V-photon corre¬ 
lation interference events. 

Even more interestingly, the same simple conditions 
lead to the computational hardness of the exact MBCS 
problem, establishing a connection between the occur¬ 
rence of multiphoton correlation interference and com¬ 
plexity. Indeed, for approximately equal detection times 
td~tGT and equal polarizations pd = e\, ^d G T), the 
multiphoton detection probabilities in Eq. ( 6 ) become 


p(V,S) 

^{td,Pd} 




sGS 


The interference of all A^-photon quantum paths in 
Eq. (9) depends, apart from an overall factor, only on 
the permanent of a submatrix of the interferom¬ 

eter random unitary matrix U. For N ^ M, these ma¬ 
trices have elements given by approximately independent 
and identically distributed (i.i.d.) Gaussian random vari¬ 
ables and the approximation of their respective perma¬ 
nents is a ^P-hard task [14]. We emphasize that the 
presence of only an arbitrarily small fraction of samples 
with probabilities as in Eq. (9) would be enough to en¬ 
sure the hardness of the exact MBCS. This can be shown 
analogously to the hardness proof of the original prob¬ 
lem of exact boson sampling in [14] . Indeed, the ability 
to perform exact MBCS with a polynomial number of 
resources would imply that the task of approximating 
any given, fixed permanent associated with the proba¬ 
bility distribution ( 6 ) is in the complexity class BPP^^. 
Since this would also include the task of approximating 
the ^P-hard permanents emerging in Eq. (9), the poly¬ 
nomial hierarchy would collapse to the third level, which 
is strongly believed to be highly unlikely. We refer to 
section II of the Supplemental Material for more details 
[57] . 

Interestingly, differently from the original boson sam¬ 
pling problem [14], the classical intractability of exact 
MBCS is not conditioned on input photons with approx¬ 
imately identical spectra in Eq. (2). Only the simple 
conditions (5) and ( 8 ) on the spectra are enough to guar¬ 
antee its computational hardness. 


Approximate MBCS. Is approximate MBCS also 
not tractable with a classical computer? Such a ques¬ 
tion is obviously of fundamental importance from an ex¬ 
perimental point of view, since it takes into account the 
inevitable experimental errors in an MBCS quantum in¬ 
terferometer which make only approximate sampling pos¬ 
sible [58] . We consider, for simplicity, the case of an N- 
photon interference matrix in Eq. (7) with unit elements 

a(s,s') = l Wsjs'gS. ( 10 ) 


This corresponds to two possible scenarios. Either all the 
input photons are completely identical or they differ only 
by their color, i.e. central frequency. In these cases the 
input photons have equal polarizations and are always 
indistinguishable at the detectors independently of the 
detection times and polarizations. 

To simplify the expressions, we consider here 
polarization-insensitive-detectors. 

a. Identical input photons. For approximately iden¬ 
tical frequency spectra 

^s{uj)=^{uj) Vs e 5, 


by using Eq. ( 6 ), the polarization-insensitive detection 
probability reads 


P, 


(V,S) 


{td} 


{td^Pd} 

{Pd}^{ei,e2}^ 

peTmU<-^’^^\\2Tif [] \xitd)\\ 

d^V 


where we used the property I]prf=ei.e 2 \Pd ' ^ 

Of course the only possible events occur within a 
detection-time interval where the function |x(t(i)| = 
l ■^[^](0 — V\t)| is not negligible. Here, independently 
of the detection times {tdjdg-p, all the probability rates 
associated with each possible sample {T>^{td}de:v) are 
given, apart from a prefactor, by the permanents oiNxN 
submatrices of the interferometer transformation 

U. 

When the observer ignores the information about the 
detection times the approximate MBCS problem reduces 
to the well known standard formulation of the approx¬ 
imate boson sampling problem, which Aaronson and 
Arkhipov argued to be intractable with a classical com¬ 
puter [14] . Therefore, the approximate MBCS problem 
is at least as complex as the original approximate boson 
sampling problem. 

b. Photons of different colors. We now address the 
case of input photons in Eq. (1) with spectral distribu¬ 
tions 


^s{oj) = u^(a;-a;s)e'“‘«, 

with equal emission times tos = and equal polariza¬ 
tions Vg — V but different colors ujg. For simplicity, we 
consider spectral shapes 
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with equal bandwidths Aujs = Auj < jo’s—Ws'| Vs, s', 
where sine a: — sina:/a;. The 7V-photon interference at 
the detectors is therefore characterized by the Fourier 
transforms 




Xs{t) = vd — lect ( 


with the rectangular function 


1 N<2 

rect a; := < 1 i i 1 

2 I""! = 2 

0 else 


Therefore, the condition (10) is satisfied, and the proba¬ 
bility rates in Eq. (4) are non-vanishing only for detection 
times td&T — [tQ + At—l/Auj, 1^ + At+1/Auj] Vd G V. 
Moreover, Eq. (5) is fulfilled for integration times 


Tf «C Iws — Ws'I ^Vs,s'g5, (11) 


where {TiAuj)~^ defines the number of discrete steps 
of length 2Tj along the time interval T. As is known, 
detectors with such high time-resolution cannot distin¬ 
guish photons of different colors Wg and multiphoton in¬ 
terference can be observed. Indeed, from Eq. (6), the 
polarization-independent detection probabilities are 


Pllf = (A^T,) 


N 


perm I [U^ 


([' 




d.s 


dGV 

sGS 


( 12 ) 


for all possible detection time intervals [td — T],td + T/] C 
T. Such probabilities are proportional to permanents 
of matrices whose elements are the elements of 
multiplied by the complex phases expiiujstd)- 

Since the elements of the submatrices are 

i.i.d. Gaussian random variables and the phase factors 
only rotate such elements in the complex plane, the 

entries of the matrices are also i.i.d. 

s^S 

Gaussian random variables as shown in App. G of the 
Supplemental Material [57]. Therefore, the probability 
distribution of the interferometer output interestingly de¬ 
pends, for all possible samples, on permanents whose ap¬ 
proximation to within a multiplicative factor is a ^P- 
hard problem [14]. Gonsequently, even for input photons 
of different colors, it is possible to show in analogy with 
Ref. [14] that approximate MBGS is of at least the same 
complexity as the standard boson sampling with identi¬ 
cal photons [60] . As a “bonus”, the number of possible 
samples (V, {td}de-D) is exponentially larger (by a factor 
{TiAuj)~^ with TjAuj ^ 1 according to Eq. (11)) with 
respect to the standard boson sampling problem. 

Does approximate MBGS retain its complexity even for 
photons which are completely pairwise distinguishable in 
their colors Wg (i.e. Wg — Wg/ Aoj Vs ^ s')l We first 


emphasize that, since these photons are characterized by 
a pairwise overlap 

OO 

J duj^si^^) ■ = 0 (13) 

0 

the approximate boson sampling problem is trivial [14] . 
Indeed, in this case, by averaging the rate in Eq. (4) over 
all possible detection times and polarizations, one finds 
that the boson sampling probability [1] 

= perm |^]dec, 

’ sGS 

for an output port sample V, is given by the permanent of 
a non-negative matrix that can be approximated with a 
polynomial number of resources [61] . Gonsequently, one 
might guess that also the approximate MBGS is compu¬ 
tationally trivial. Nonetheless, the complexity emerging 
from the result in Eq. (12) is independent of the colors 
Wg of the input photons, demonstrating that also in this 
case approximate MBGS is classically intractable. 

Two essential physical aspects are behind the demon¬ 
strated complexity of approximate MBGS: all possible 
detection-time events can be an outcome of the sam¬ 
pling experiment (none of the events is disregarded) and 
all these time samples arise from the interference of A^! 
multiphoton quantum paths. In conclusion, the physics 
of sampling among all possible Wphoton interference 
events behind our proposal is at the heart of the com¬ 
plexity of approximate MBGS. 

Discussion. In this letter, we demonstrated how and 
to what degree the occurrence of multiphoton interfer¬ 
ence in time- and polarization-resolving correlation mea¬ 
surements leads to computational hardness in linear op¬ 
tical interferometers. 

The definition of an A^-photon interference matrix 
a{s, s') in Eq. (7) allowed us to formulate the simple suffi¬ 
cient condition (8) on the spectra of the input photons for 
the occurrence of A-photon interference, provided suffi¬ 
ciently small integration times (see Eq. (5)). 

Remarkably, these two simple conditions are also suf¬ 
ficient to guarantee the complexity of exact MBGS. In 
contrast, the complexity of the original exact boson sam¬ 
pling problem has only been proven for identical input 
photons. 

For approximate MBGS on the other hand, not only 
the existence of samples exhibiting full A-photon inter¬ 
ference (guaranteed by (8)) is important but also their 
fraction with respect to the total number of samples. In¬ 
terestingly, this is encoded in the magnitude of the entries 
a{s, s') of the A-photon interference matrix in Eq. (7). 

It was thus natural to consider the simple case of full 
overlap of the modulus of the single-photon detection am¬ 
plitudes (a(s, s') = 1) where all possible detection events 
correspond to A-photon interference samples [62] . In this 
case, corresponding to identical input photons or photons 
with arbitrary colors, approximate MBGS is at least of 
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the same complexity as boson sampling with identical 
photons. 

This is particularly interesting if the differences in the 
central frequencies are much larger than the width of the 
single photons’ spectral shapes, corresponding to fully 
distinguishable photons in the sense of Eq. (13). While 
approximate boson sampling becomes trivial in this case 
[14], approximate MBCS is at least as complex as when 
perfectly identical input photons are used. 

Since detectors with high temporal resolution (< Ins) 
and single photons with large coherence times (> Ips) 
are readily available today experimentally [36], the re¬ 
quirement of time-resolved measurements in the imple¬ 
mentation of the MBCS problem can be readily ful¬ 
filled. Moreover, an implementation of MBCS has the 
advantage to ease the difficulties faced in the produc¬ 
tion of identical photons. Indeed, photons of approxi¬ 
mately equal colors {Au! 3> jwg — lUs' \ Vs ^ s') are not 
needed any more, unlike in the original approximate bo¬ 
son sampling problem. This furthermore paves the way 
towards the use of photons of arbitrarily small bandwidth 
Auj, where the indistinguishability in the emission times 
(1/Z\a; \tos — ios'l Vs, s') can be easily achieved. 

In conclusion, all these results represent an important 
stepping-stone towards a full fundamental understanding 
of the complexity of multiphoton interference of photons 


of arbitrary spectra in linear optical networks, when the 
information about detection times and polarizations is 
not ignored. This may lead to “real world” applications 
in quantum information processing [52] and in quantum 
optics overcoming the experimental challenge in the pro¬ 
duction of identical bosons. 

Finally, our results can be extended to bosonic inter¬ 
ferometric networks with atoms [12, 13, 63], plasmons 
[64] or mesoscopic many-body systems [65] and are also 
relevant to the study of the complexity of multiboson cor¬ 
relation interference for different input states [2, 66, 67] 
and different correlation measurements [68] . 
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From the Physics to the Computational Complexity of Multiboson Correlation 

Interference: Supplemental Material 


I 


I. THE MULTIBOSON CORRELATION SAMPLING (MBCS) PROBLEM 


Here, we formally define the MultiBoson Correlation Sampling (MBCS) problem in linear interferometers described 
by a unitary M x M matrix chosen randomly according to the Haar measure. 

We first introduce the notion of families of interferometer input states {|1 [^s])s}sg5i at the input ports s G S, 
defined by the sets {^s}sg 5 of N complex spectral distributions 

■■= Vs^s{uJ - (SI) 

with square integrable spectral shapes — ^s) centered around uj = lOs, polarizations Vg, central frequencies Ws, 
and emission times t^s- 

Definition 1 (Families of input states). Let I be the set of square normalized, complex spectral distributions which 
fulfill the narrow-bandwidth approximation. For given sets F/v C J^, with S IN, of A-tuples {^s}sg 5 G Cv of 
single-photon spectra, we define the family F ■= {FatItvgin- Further, we denote as F the set of all possible families 

F. 


Secondly, we formally define all the possible samples that can be detected at the interferometer output in the 
MBCS problem by discretizing the detection-time axis into bins of width 2Tj (integration time) much smaller than 
the temporal widths of the photons and centered at times tk ■= 2Tik. Here, k = A:inin + where 

^min, ^max G ^ define the temporal interval where detections can occur. From now on, we will refer to the time 
samples {tdjdG'D only in terms of the N integers {kd = td/{2Ti)'\d^'D and define each possible overall sample as 
© := {V; {kd}, {Pd}), with pd G {ei, 62 } for a fixed orthonormal polarization basis {ei, 62 }. 

The probability distribution 

Vf ='PfH^s} G Fn, A G iiM,N) (S 2 ) 

associated with all the possible samples 6 is defined by a given rectangular M x N submatrix 


A ■■= [hld,e]d=l,...,M G Am,N 
S^S 

of the interferometer matrix U, where A.m,n is the set of M x TV matrices with orthogonal columns, and by a given 
set of complex spectral distributions {^s} G Fm for the N input photons, with F = {FArjArgiN G F. 

As shown in the main letter, defining the matrices 


^ {kd,Pd} 


n-FF) 

{tk^,Pd} 


l^d,s{pd ■ Xsi^Tjkd)) 


de-D’ 

sGS 


(S3) 


where 


Xs{t) = VgXsit - tos - At) e 


iujs (t—tQs — At) 


(S4) 


are the Fourier transforms of the spectra ^s(w) in Eq. (SI) (At is the delay the photons pick up in the interferometer), 
the probability pe associated with the sample 6 , when sampling from the distribution Vf Ai Eq. (S2), is 


Pe 


■■= Pfl®] = <fL, = mr permr/P'Si 


^ {tkd^PA ~ '{kd,Pd} 

For samples ©id with identical detection times kd = k Wd G V and polarizations pd = p G V, this simplifies to 


(S5) 


PSm = i‘2Ti)^ n \P ■ Xs(2Tik)f IpermT^^®’'^) 

sGS 


(S 6 ) 


Following Aaronson and Arkhipov in [SI], it is reasonable to represent the elements of the interferometer matrix as 
rational numbers (x+iy )with integers x and y. Analogously, the same holds for the values T/, Xs(td—tos—At), 
Pd ■ Vg, and exp(iujg{td — tos ~ At)), with s G S (notice that this representation is reasonable since all these given 
complex values describing an MBCS experiment can be approximated with high enough precision by choosing poly(A) 
sufficiently large). 

We emphasize that, while we will prove the hardness of the MBCS problem assuming the rational representation of 
these numbers, this proof can easily be generalized to a representation in terms of algebraic numbers which are dense 
in C (more information can be found in App. B). 
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Theorem 1 (Exponential size of the probabilities). Assume that the interferometer matrix and the temporal dis¬ 
tributions are given in a rational representation, as described above. Then, independently of the specific choice of 
A G Um.n! F = {E/v} G F, and {^s}sg5 € Fjq, all the probabilities (S5) in the probability distribution Vp o,t most 
exponentially small, 

pe > V6. 

Alternatively, the same result holds if we instead assume that the values are represented by algebraic numbers (which 
are dense in G). 

Proof. This is demonstrated in App. A. ■ 

Furthermore, the probabilities can be represented as (see App. A) 

_ We 
P® “ 2P°h(w) 

with an integer we, which ensures that an MBCS oracle as it will be defined in Definition 2 (using a random input 
string r of length at most polynomial in N) is able to sample from a probability distribution equal or arbitrarily close 
to the exact probability distribution in Eq. (S2) where the probability for each sample 6 is defined by Eq. (S5). 

We can finally define an MBCS oracle as: 

Definition 2 (Definition of an MBCS oracle in analogy to [SI]). Let Op be an oracle that takes as input an M x N 
matrix A € Hm.Nj an Wtuple S Fn from a given family F G F, an error bound /3 > 0 (encoded in 

unary as 0^^^ to ensure a scaling of resources with 1//3)), and a string r G {0,1}^°^^^'^^ which is its only source 
of randomness. Let = P’of({€«}) A;/3) be the distribution over the outputs of Op if A, {|s}) and f) are 

fixed but r is uniformly random. Then, Op is called an exact MBCS oracle for the family F if 'PofU^s}, A; /3) = 
7^f({€s}, A) for all G IN, A G Hm.tv, {^s} S Fn, and j5 > 0. Further, Op is an approximate MBCS oracle for 
the family F if \\Vof{{^s}, A; jd) — 'P_f({^s}, A)|| < fd for all G IN, A G Hmn, {Is} G Fn, and (d > 0, where 
\\VoFim,A-,p)-Vpi{^s},A)\\ ■■=i/2j:^\PTp^je]-PTvAe]\- 

II. THE COMPLEXITY OF THE EXACT MBCS PROBLEM 

Here, we formally demonstrate the complexity of the exact MBCS problem for the set £ of families E defined as: 

Definition 3 (Set £ of “complex” families in the exact MBCS problem). We define £ C F (with F defined in 
Definition 1) as the subset of families E ■■= {En}n^n which fulfill the condition 

VAf G ]N,V{|s} G Eat : 0<a(s,s')<l Vs, s' G 5, (S7) 


with a(s,s') := \vs ■ Vs'\J(^^dt |xs(t- tos)| |Xs'(t - tos')l- 
We first introduce the following theorem: 

Theorem 2. Let E = {En}n€:N & £■ It is then ensured that a sequence of time-polarization tuples, 

with Un G {kmim • • ■) kmax\ ond pn G jei, 62}, exists such that 

\Pn ■ Xs{2TikN)f > 0 Vs G 5. 

For any given N, the tuple {kN,PN) can be found in polynomial time. 

Proof of Theorem 2. As already discussed in the main letter, the condition (S7) for a fixed N directly implies the 
existence of at least one time bin /cat and one polarization pn for which all N amplitudes pn • XsiflTikN), with s G S, 
are non-vanishing. 

Moreover, due to the finite number L ■= fc^ax ~ ^min + 1 of time bins in the interval [femin, ^max]) {kN^PN) can be 
found in polynomial time. I 

We can now formulate the main theorem on the complexity of exact MBCS: 

Theorem 3 (Main theorem on the complexity of exact MBCS). If Op is an exact MBCS oracle for a given family 
E G £, then C BPP^^ . This implies that the polynomial hierarchy collapses to the third level if the exact 
MBCS problem for states in the family E can be solved in polynomial time by a classical computer. 
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Proof of Theorem 3. If a matrix X G and a parameter g € [1 + l/poly(-/V), poly(7V)] are given, it is ^^P-hard 

to approximate |permX| to within a multiplicative factor g [SI]. 

We now show that, given an MBCS oracle Oe for a family E G £ (with £ defined in Definition 3), it is possible to 
perform this approximation in FBPP’^^ ^. 

As shown in [SI], given 7 := l/|jW|| > it is possible for any M > 2N to find in polynomial time a unitary 

M X M matrix U which contains 7 X as its top-left N x N submatrix, i.e. 

'yX = [Ud,s]d=i,...,N- 
s=l,...,N 

Given S — {1, 2,..., N}, the matrix A = [Z^d,s]d=i,....M G Am,n and the spectra {4 s}sg5 G E^r induce the probabil- 

s^S 

ity distribution Ve = 'Pe{{^s}, A)- From Theorem 2, a tuple (fc 7 v,PAr) exists such that JIsgS \Pn • Xs{2TjkN)\ > 0. 
For the sample 


6* := {V*]{kN, • ■ ■, k^}, {pn, ■ ■ ■,Pn}), 


N times 


N times 


with 


V* :={1,2,...,7V}, 


the probability in Eq. (S 6 ) becomes 

PS* = Pr[6*] = n \pN ■ Xs{tN)\" jpermAj' . (S8) 

® sG 5 

According to Definition 2, the exact MBCS oracle Oe samples from the probability distribution VoeH^s}, A-, P) = 
P£;({^s}, A) returning the value 0£;({|s}, A,r) = 6 * with probability 

pe*= ^ Pi [OEim, At) = e*]. 


Defining the Boolean function 




OE{m,A,r) = e* 
OE{m,A,r)y^e* ’ 


we can also express this probability as 


Pe* 


1 

2 poly(Ar) 


n /W- 

rG{0,l}P°'y<”) 


(S9) 


(SIO) 


Therefore, we can use Stockmeyer’s algorithm [S2[ to approximate ps* to within a multiplicative factor ps[l-l- 
l/poly(A^),poly(A^)] in FBPP'^p ^ in time polynomial in N. Consequently, from Eq. (S 8 ), we can approximate 
jpermAj^ in FBPP^p ^ as well. Since performing such an approximation is a ^(tP-hard problem [SI], P#^ C 
gppNP E MBCS problem could be solved in polynomial time by a classical computer this would imply 

P"^^ C BPP^^, which, by Toda’s theorem [S3], would lead to a collapse of the polynomial hierarchy to the third 
level. ■ 


III. THE APPROXIMATE MBCS PROBLEM 

We consider, for simplicity, a multiboson correlation experiment with polarization-insensitive detectors and states 
from a family R := {RN}NeN G E (with E defined in Definition 1) of input states {|l[^s])s}sg 5 defined by input 
spectra 


{4s(w)}sgS 


/ 1 . 

< V , sine 

t V ttAuj 


LU — LJs 

Aid 


^iuto 


sGS 
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with equal emission times tos = to, equal polarizations Vg = v, and equal bandwidths Aujg = Auj 
but different colors ojs. For these spectra, the elements of the A^-photon interference matrix a(s,s') ■■= 
ks • ■Ws'l dt \xs{t- tos)| \Xs'{t - tos')l fulfill the condition 


a{s, s') ~ 1 Vs, s' 

of full temporal overlap. As in section II, we discretize the detection-time axis into time bins of width 2Tj. Here, 
as described in the main letter, in each output port a detection can only occur in the time bins tk ■■= 2Tik G T ■■= 
[to + At - l/ALo,to + At + l/Auj], with k e fcmin + 1, • ■ •, fcmax}, where k^in ■= {to + At - IjAw + Ti)/{2Ti) 

and fcniax •= (to + At -I- 1/Alo — Ti)/{2Tj). Therefore, the total number of time bins in which a detection is possible 
is L ■■= l/{TjAui). Given the output probability distribution 


'Pr ■= ’^ii({^s}sG5 V Rn, a e ilM.w), 

for a given sample © — (P; {fed}), the probability 


pe := Pr[6] = 4^5) = {TjAuj)^ perm 


derived in the main letter, takes the form 


Pe 


perm 



^2iTiu}^kd 



(Sll) 


We can now formulate the following theorem on the computational power of an approximate MBCS oracle for the 
family R\ 

Theorem 4 (Computational power of an approximate MBCS oracle for the family R of input states). Let Or be an 
approximate oracle for the MBCS problem for the family R and e,5 > 0 given error bounds. It is then possible to 
approximate the modulus square of the permanent of a random Gaussian N x N matrix to within an additive error 
zteNl and with success probability larger than 1 — 6 in FBPP*^^ ^ in time polynomial in N, Ije, and 1/(5. 

Proof of Theorem /. Let A be a complex N x N matrix whose entries are randomly picked according to a standard 
complex normal probability distribution ^(0, l)c. We will adapt the arguments found in [SI] to give an algorithm in 
FBPP^^ ” that performs the required approximation in polynomial time. The main idea is to introduce axi M x N 

matrix A = [Ud,s]d-i,...,M G Am,n (as in the exact case, A is a rectangular submatrix of a Haar-random unitary 
s^S 

interferometer matrix IL corresponding to the input configuration S ■■= {1,2, ...,A}) such that a specific sample 
©* := {'D*,{k’^}) occurs with probability pe* oc |permA| , which can be estimated using the approximate MBCS 
oracle defined in Definition 2. 


1. In order to rule out the possibility that the oracle can willingly sabotage the output probability of the sample 
©* in which we are interested, ©* has to be picked randomly. Therefore, the time sample {k^} is generated 
by randomly picking N time bins tk* = 2Tik'^ € T according to a uniform probability distribution. Given the 
structure of the probabilities in Eq. (Sll), it is useful to define a matrix X which incorporates the inverse of 
the phase factors corresponding to the randomly chosen time sample {k^}, i.e. 




— 2iTi(jJs - fc j 


where di and Si are the f-th elements of T) and S, respectively. Since the elements of X are i.i.d. Gaussian 
random variables, the same holds for the elements of X, as shown in App. C. 

2. An M X N matrix A € Am.n is generated with the following properties: A is distributed like the submatrix 
{fdd,s]d=i,...,M of a Haar-distributed M x M unitary matrix U and it contains X' ■= X/'/M as a uniformly 

random submatrix. As shown in [SI], since A is a Gaussian matrix this is possible to achieve in BPP'^^ with a 
failure probability 


Pr [hiding failed] < 


6 

4’ 


(S12) 






V 


provided 


N 

M = K— log^ — (S13) 

0 0 

with a sufficiently large constant K. The random output-port sample that corresponds to the position of 
X' = [%,s]dG’D* inside of A is denoted as T>*. Thus, for the random sample ©* — {'D*\{k*^}), Eq. (Sll) 

S^S 

becomes 


Pe* 


1 


perm 



3 = 1,-X. 


2 


^l4^|permX|Z 


(S14) 


3. While pe* describes the probability for the sample ©* within the exact probability distribution Vr = 
'Pr{{^s},A), we define 


as the probability for the sample ©* within the probability distribution Vor = 'Pouii^s}, A; f3) of the approx¬ 
imate MBCS oracle On (as shown in App. A, all probabilities pq are at most exponentially small and it is 
thus ensured that an oracle On as defined in Definition 2 can sample according to a probability distribution 
^Oh({€s}, A;/3) arbitrarily close to Vnil^s}, A)). 

In analogy to Eqs. (S9) and (SIO), we can also define a Boolean function 




On{m,A,r;(3) = e* 

On{{is},A,r;(3)^e*^ 


and write 


ge* = 


1 

2 poly(Ar) 


H /W- 


This makes it evident that Stockmeyer’s algorithm [S2] can be used to find a value gs* in FBPP^^°” that 
approximates ge* to within a multiplicative factor (1 -I- a) S [1-1- l/poly(iV), poly(-/V)], such that 


Pr [ge*/(l + «) < ge* < (1 + a)ge*] > 1 - ^ 


(S15) 


in time polynomial in M and 1/a [S2]. 


4. By using Definition 2 of the approximate MBCS oracle and Eq. (S15) for a = £(5/16, we show in App. D that 
the estimate qe* of the squared permanent of X satisfies the condition 


Pr 


IpermAr 


> eN\ 


6 

^ 2 


1 


(S16) 


Adding to Eq. (S16) the probability of failure for the hiding of X' from Eq. (S12) and recalling Eq. (S13) which 
ensures 2~^ < (5/4, we find that the total probability of failure is smaller than 6, as required. 

Further, the algorithm runs in a time polynomial in N, 1/6, and 1/e since the hiding procedure in step 2 is 
polynomial in time with respect to N and 1/5 and the running time of Stockmeyer’s algorithm (used in step 3) is 
polynomial in N and 1/a, where a = ei5/16 was chosen. ■ 

In [SI], the authors argue that approximating the modulus square of the permanent of a Gaussian matrix is a 
/^P-hard problem if two reasonable conjectures are true. Under this assumption, the following theorem holds: 

Theorem 5. Let On be an approximate MBCS oracle for the family R. If the approximate MBCS problem for R can 
be solved in polynomial time by a classical computer, the polynomial hierarchy collapses to the third level. 

Proof of Theorem 5. Theorem 4 states that it is possible to approximate the modulus square of a permanent in 
PBPpNP given an oracle On for the family R. If this approximation is indeed #P-hard, as argued in [Sl[, it 
follows that P#^ C BPP^^ . Further, if the MBCS problem for states of this family can be solved in polynomial 
time by a classical computer, then P^^ C BPP^^ which implies, by Toda’s theorem [S3] , that the polynomial hierarchy 
collapses to the third level. ■ 
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Appendix A: Exponential lower bound on the MBCS probabilities 

Here, we show that the assumption of a rational representation (a: + *j/)/ 2 P°h(-^v) integers x, y) of the elements 

of the interferometer matrix and of the values of the temporal distributions for all possible time bins results in an 
exponential lower bound on the non-vanishing detection probabilities pq . 

(T) 

Using the expression (S5), the definitions of a matrix permanent and of the matrix in (S3), and the 

explicit expression Eq. (S4) for the functions Xs(t)i we find 

pe = \ n ■ ^.'{d^TiPd ■ 

(7,a' ^ d^X> 

X n (Xa'id){td - toa'id) - ^t))*Xa{d){td - tocr(d) “ 
den 

X (d)('td — to„Hd)— Qi'^a{d)('td — toa{d)—‘^t) 

dec 

Since all 8A -|- 1 numbers Tj, lAd,s and c.c., {pd ■ v^) and c.c., Xsi^d — tos — and c.c., exp[ia;s(td — tos — 2\t)] and 
c.c. are represented as (x -|- z?/)/ 2 P°h(w) integers), it is immediately clear that the probability has the form 

1 

“ 2poly(7V) 

with a non-negative integer we and is therefore at most exponentially small. This guarantees that Stockmeyer’s 
algorithm can be used to approximate these probabilities. 

The probabilities are also at most exponentially small if these values are represented by algebraic numbers. Indeed, 
it was shown in [S5] that 


Pe 


> 


2-r(N) ^ 


(SAl) 


with a polynomial r{N). 


Appendix B: Using algebraic numbers instead of rational numbers with polynomial precision 

Here, we will show that, instead of assuming that the values lAd,s, Ti, Xsitd — tos — ^t), and exp[iujs{td — ~ ^i)] 

are represented as rational numbers {x + iy )as in [SI], the hardness of MBCS sampling can also be proven if 
these values are represented by algebraic numbers [S4] . This is appealing because it is possible to find an arbitrarily 
close algebraic approximation of any complex number since the algebraic numbers are dense in C. 

We emphasize that, in general, an MBCS oracle as defined in Definition 2 now cannot sample from the exact 
probability distribution Vf in Eq. (S2) any more. However, as we will demonstrate in the following, the hardness 
proof in section H is still valid if we define an “exact” MBCS oracle as an oracle which samples from a probability 
distribution Vop with 


ro.({4.},A)-7^^({6},A)|l<2-W, 


(SBl) 


where the polynomial s{N) is assumed to dominate the polynomial r{N) from Eq. (SAl) {s{N) — r{N) grows mono- 
tonically, and s(l) > r(l) -|- 2). Such an approximation of the exact probability distribution can however still be 
achieved by an oracle as defined in Definition 2. 

Defining the probabilities qe — Pr^^^]©], Eq. (SBl) implies that 


be-gel <2-2-*W 

PS - 2 • < ge < Pe + 2 • 2-"W 


1 - 2 - 


2-s{n) 

Pe 


< ge < Pe 1 + 2 


2-s(n) 

Pe 


^ Pe 
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which with Ineq. (SAl) becomes 


Pe 


Pe 1 - 2 


2-s(n) 

2-Hn) 


^ 96 ^ Pe 1 + 2 


2-s(7V) 

’ 2 -+^) 


(l _ 2-+W-+^)-i)) < ge < P 6 (1 + . 


Using the inequality 1 — x > 1/(1 + 2x) for 0 < x < 1/2, we find that 

1 


Pe 


1 + 2 -+W-’’W- 2 ) 


<96<P6 (i + 2-(*W-’'W-2)), 


(SB2) 


i.e. that qq is multiplicatively close to pQ to within a factor 1 + 2 

The approximation of pq to within a multiplicative factor g G [1 + can therefore be achieved 

by approximating qq to within a factor g' ■■= g/(l + Since g — 1 is at most polynomially small, 

g' G [1+ poiy(jv) ’ PO^y(-^)] well. Therefore, Stockmeyer’s algorithm tells us that a value qe can be found in BPP^^, 
such that 


9 ' 


<qe< g'qe ■ 


Then, Ineq. (SB2) implies that 


Pe 

9 


Pe 


g'{l + 2-(+^)-U1v)-2) 


<96<P65'(1 + 2-(*W-’'W-^)) 


Pe9- 


and the arguments in the hardness proof for the exact MBCS in section II still hold if we define “exact” in the sense 
of Ineq. (SBl). 

Further, the proof for the hardness of the approximate MBCS problem from section III is still completely valid in the 
algebraic number representation since an approximate MBCS oracle as defined in Definition 2 is still able to sample 
from a probability distribution that is polynomially close in variation distance to the exact probability distribution. 


Appendix C: The distribution of the elements of X 


We will now show that, under the assumption that the elements Xij, i, j = 1,..., iV, of a complex N x N matrix 
X are i.i.d. random variables with a complex standard normal distribution, the same holds for the elements of X 
which are defined as 


v. . — X p+'.j 

If the elements of X are i.i.d. N{0, l)c variables, their joint probability distribution is 

^ 1 

/x({V,,,})= l[ 

j. j. 7]- 

* j'=l 

Noting that the Jacobi determinant for the change of variables between the {Xij} and the {Xij} is det J = 1, we 
find that the common probability distribution for the elements of X is 

N ^ 

/+(fe}) = /x({A,,7e-‘^^..}).^tJ= n 

—V— TT 

= 1 *,i = i 


where we used that the complex Gaussian distributions are independent of the complex phase. 

Therefore, the elements of the matrix X are still i.i.d. ^(0, l)c random variables if the elements of X were. 
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Appendix D: Bound on failure probability of approximation 

As in [SI], we define 'Pm,n as the set of all port samples V and Gm,n as the set of bunching-free port samples, i.e. 
as the set of port samples T> which consist only of pairwisely different output port indices. Further, let be the 

set of all time samples {kj}- 

We will now proceed to derive three inequalities which combined yield Ineq. (S16). 


1. With Aq ■■= \pq — gel, we find the expectation value 


aI 


< 


def AI'POr —'PrW ‘^W'POr —'Pr\ 


\Gm,n \ 


^ (")i" 


\Gm,n\ \^l,n\ 

<3^ 


(LM) 


N ’ 


where, in the penultimate step, we used the Definition 2 of an approximate MBCS oracle and, in the last step, 
the fact that M = uj{N'^) {M asymptotically grows faster than N'^) [SI]. Therefore, Markov’s inequality gives 
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SGGm.niSi W, 


A© > 3/3 


Nl 


{LMY 5 


5 

<4 


which, by choosing /3 = £(5/24, simplifies to 


Pr 

&€Gm,N'S’^L,N 


A© > 


£ m 

2 (LM)^ 


< 


(SDl) 


The MBCS oracle only knows A but has no information whatsoever about which sample 6* has been chosen 
from Gm,n <8 Thus, Eq. (SDl) implies 


Pr 

X,A 


A©. > - 


N! 


2 (LM)^ 


S 

<4- 


(SD2) 


2. As shown in step 3 of the proof of Theorem 4, we can apply Stockmeyer’s algorithm [S2[ to find an approximate 
g©* of g©* in FBPP^^ The algorithm guarantees that, for an arbitrary a > 0 and a runtime polynomial in 
1/a and M, 


Pr[|g©. - g©.| > ag©.] < Pr[g©. > (1 + a)g©. V g©. < g©*/(l + a)] < 


2M’ 


where the first inequality follows from 1 — a< 1/(1 +a) and the second inequality is equivalent to Ineq. (S15). 
With the choice a = £i5/16, this becomes 


Pr[|g©, -g©.| > g©.] < 


(SD3) 


3. Lastly, 


f 1 — 0^^L.N ^*5 ^ _1_ _ 

SGGM,N®Gi,,iv |Gm,7v| \Gl,n\ ~ IGm.AtI \^L,n\ Cn)L^ (LM)^ 


1 fV! 

< 2 - 


and thus by invoking Markov’s inequality 
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(LM)2 <5 

With the same arguments leading to Ineq. (SD2), this implies 
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IX 


Indeed, by using the inequalities 

Pr[|a — 5| > c] < Pr 

and 

Pr[|a — 5| > c] < Pr 
demonstrated in App. E, we obtain 
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I I c 

|a| > 2 


Pr 


\b\>^\ (c>0) 


^>kl 


+ Pr [|a — 5| > kb] (c, k > 0), 


(SD5) 


(SD6) 
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2 (LM)^ 


4 2^ 4 2 2^' 


where in the last step, we inserted Ineqs. (SD2), (SD3), and (SD4). This inequality is equivalent to Ineq. (S16), as 
can be seen by inserting the expression for from Eq. (S14). 


Appendix E: Inequalities for probabilities 


Here, we want to prove the inequalities (SD5) and (SD6). 
Proof. First, observe that 


and therefore 


Pr 


|a|<|A|5|<|^|a-&|<c 


|a|<^A|6|<^l <Pr[|a-&|<c]. 


2 ' ' - 21 

Second, for two events A and B, we know that 

Pr[A nB]= Pt[A] + Pr[B] - Pr[A U B] 
and it follows that {A, B denote the complement of A, B) 

1 - Pr[A n H] = 1 - Pr[A] - Pr[H] + Pr[A U H] < 2 - Pr[A] - Pr[H] = Pr[A] + Pr[i?]. 
Now, Ineq. (SD5) follows as 


Pr[|a — 6| > c] = 1 — Pr[|a — 6| < c] < 1 — Pr |a| < ^ A |6| < ^ 


< Pr 


I I c 

l“l > 21 


Pr 


1^1 > i; 


where we used Ineq. (SEl) and Ineq. (SE2) in the second and third step, respectively. 
Ineq. (SD6) can be proved in a similar way. 


(SEl) 


(SE2) 
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